function [snapshot]=afd_bc_outer(delx,delt,velocity,snap1,snap2,snapshot,boundary)
% AFD_BC_OUTER ... apply absorbing boundary condition to outer boundary
%
% [snapshot]=afd_bc_outer(delx,delt,velocity,snap1,snap2,snapshot,boundary)
% AFD_BC_OUTER applies absorbing boundary conditions to the outer layer of rows
% and columns in a situation where there are 1 or 2 layers of rows and columns on
% the boundary.  The output matrix has absorbing boundary conditions along its
% outer edge. In order to be stable, vmax*delt/delx must be < 0.7.  
%
% delx = the grid spacing for horizontal and vertical bins in consistent units
% delt = the time interval in seconds
% velocity = the input velocity matrix in consisnent units
%          = has a size of floor(zmax/delx)+1 by floor(xmax/delx)+1
% snap1 = the wavefield at time=t-2*delt (same size as velocity matrix)
%        = will be based on the source array desired i.e. the position
%          of the sources will be one, and the rest of the positions
%          will be zero
% snap2 = the wavefield at time = t-delt(same size as velocity matrix)
% snapshot = the wavefield at time = t (same size as velocity matrix)
% boundary = indicate whether all sides of the matrix are absorbing
%          = 0 indicates that no absorbing boundaries are desired
%          = 1 indicates all four sides are absorbing
%          = 2 choses three sides to be absorbing, and the top one not to be
%             this enables sources to be put on the surface
%
% snapshot = the wavefield at time = t with absorbing boundary conditions
%
% by Carrie Youzwishen, Febraury 1999 
%
% NOTE: This SOFTWARE may be used by any individual or corporation for any purpose
% with the exception of re-selling or re-distributing the SOFTWARE.
% By using this software, you are agreeing to the terms detailed in this software's
% Matlab source file.

% BEGIN TERMS OF USE LICENSE
%
% This SOFTWARE is maintained by the CREWES Project at the Department
% of Geology and Geophysics of the University of Calgary, Calgary,
% Alberta, Canada.  The copyright and ownership is jointly held by
% its 'AUTHOR' (identified above) and the CREWES Project.  The CREWES
% project may be contacted via email at:  crewesinfo@crewes.org
%
% The term 'SOFTWARE' refers to the Matlab source code, translations to
% any other computer language, or object code
%
% Terms of use of this SOFTWARE
%
% 1) This SOFTWARE may be used by any individual or corporation for any purpose
%    with the exception of re-selling or re-distributing the SOFTWARE.
%
% 2) The AUTHOR and CREWES must be acknowledged in any resulting publications or
%    presentations
%
% 3) This SOFTWARE is provided "as is" with no warranty of any kind
%    either expressed or implied. CREWES makes no warranties or representation
%    as to its accuracy, completeness, or fitness for any purpose. CREWES
%    is under no obligation to provide support of any kind for this SOFTWARE.
%
% 4) CREWES periodically adds, changes, improves or updates this SOFTWARE without
%    notice. New versions will be made available at www.crewes.org .
%
% 5) Use this SOFTWARE at your own risk.
%
% END TERMS OF USE LICENSE

[nz,nx]=size(snap1);

sqrt2 = sqrt(2.0);

if boundary == 1

  % apply absorbing boundary for top 

  snapshot(1,3:nx-2) = (2.*velocity(1,3:nx-2).*delx*delt.^2)./(delx+velocity(1,3:nx-2).*delt).*...
         ((snapshot(2,3:nx-2)./(2*delt*delx) - snap1(2,3:nx-2)./(2*delx*delt) + ....
         snap1(1,3:nx-2)./(2*delt*delx)) + ...
    (-1./(2.*delt.^2.*velocity(1,3:nx-2))).*...
         (-2.*snap2(1,3:nx-2) + snap1(1,3:nx-2 ) -2.*snap2(2,3:nx-2)+ ...
         snap1(2,3:nx-2) + snapshot(2,3:nx-2)) + ...
    velocity(1,3:nx-2)./(4*delx.^2).*...
         (snapshot(2,4:nx-1) + snap1(1,4:nx-1) + snapshot(2,2:nx-3) - ...
          2.* snapshot(2,3:nx-2) - 2.*snap1(1,3:nx-2) + snap1(1,2:nx-3)));
end

% apply absorbing boundary for bottom

snapshot(nz,3:nx-2) = -2.*delx.*delt.^2.*velocity(nz,3:nx-2)./(delx+velocity(nz,3:nx-2).*delt).*...
       ((-snapshot(nz-1,3:nx-2)./(2*delt*delx) - snap1(nz,3:nx-2)./(2*delt*delx) +...
       snap1(nz-1,3:nx-2)./(2*delt*delx)) +...
  1./(2*delt.^2.*velocity(nz,3:nx-2)).*...
       (-2.*snap2(nz,3:nx-2) + snap1(nz,3:nx-2) + snapshot(nz-1,3:nx-2) - ...
       2.*snap2(nz-1,3:nx-2) + snap1(nz-1,3:nx-2)) + ...
  (-velocity(nz,3:nx-2)./(4*delx.^2)).* ...
       (snapshot(nz-1,4:nx-1) - 2.*snapshot(nz-1,3:nx-2) + snapshot(nz-1,2:nx-3) + ...
       snap1(nz,4:nx-1) - 2.*snap1(nz,3:nx-2) + snap1(nz,2:nx-3)));

% apply absorbing boundary for right hand side

snapshot(3:nz-2,nx) =  -2.*delx.*delt.^2.*velocity(3:nz-2,nx)./(delx+velocity(3:nz-2,nx).*delt).*...
       ((-snapshot(3:nz-2,nx-1)./(2*delt*delx) - snap1(3:nz-2,nx)./(2*delt*delx) + ...
       snap1(3:nz-2,nx-1)./(2*delt*delx)) + ...
  1./(2*delt.^2.*velocity(3:nz-2,nx)).*...
       (-2.*snap2(3:nz-2,nx) + snap1(3:nz-2,nx) + snapshot(3:nz-2,nx-1) - ...
       2.*snap2(3:nz-2,nx-1) + snap1(3:nz-2,nx-1)) + ...
  (-velocity(3:nz-2,nx)./(4*delx.^2)).* ...
       (snapshot(4:nz-1,nx-1) - 2.*snapshot(3:nz-2,nx-1) + snapshot(2:nz-3,nx-1) + ...
       snap1(4:nz-1,nx) - 2.*snap1(3:nz-2,nx) + snap1(2:nz-3,nx)));


% apply absorbing boundary for left hand side

snapshot(3:nz-2,1) =(2.*velocity(3:nz-2,1).*delx*delt.^2)./(delx+velocity(3:nz-2,1).*delt).*...
       ((snapshot(3:nz-2,2)./(2*delt*delx) - snap1(3:nz-2,2)./(2*delt*delx) + ...
       snap1(3:nz-2,1)./(2*delt*delx)) + ...
  (-1./(2*delt.^2.*velocity(3:nz-2,1))).* ...
       (-2.*snap2(3:nz-2,1) + snap1(3:nz-2,1) + snapshot(3:nz-2,2) - ...
       2.*snap2(3:nz-2,2) + snap1(3:nz-2,2)) + ...
  (velocity(3:nz-2,1)./(4*delx.^2)).* ...
       (snapshot(4:nz-1,2) - 2.*snapshot(3:nz-2,2) + snapshot(2:nz-3,2) +...
       snap1(4:nz-1,1) - 2.*snap1(3:nz-2,1) + snap1(2:nz-3,1)));


%%%%%% Absorbing boundary conditions for corners

% for lower right hand corner

snapshot(nz-1,nx) = velocity(nz-1,nx).*delt.*delx./(2.*velocity(nz-1,nx).*delt + sqrt2*delx).*...
       (snapshot(nz-2,nx)./delx + snapshot(nz-1,nx-1)./delx + ...
       sqrt2/(velocity(nz-1,nx).*delt)*snap2(nz-1,nx));

snapshot(nz,nx-1) = velocity(nz,nx-1).*delt.*delx./(2.*velocity(nz,nx-1).*delt + sqrt2*delx).*...
       (snapshot(nz-1,nx-1)./delx + snapshot(nz,nx-2)./delx + ...
       sqrt2/(velocity(nz,nx-1).*delt)*snap2(nz,nx-1));

snapshot(nz,nx) = velocity(nz,nx).*delt.*delx./(2.*velocity(nz,nx).*delt + sqrt2*delx).*...
       (snapshot(nz-1,nx)./delx + snapshot(nz,nx-1)./delx + ...
       sqrt2/(velocity(nz,nx).*delt)*snap2(nz,nx));

% for lower left hand corner

snapshot(nz-1,1) = velocity(nz-1,1).*delt.*delx./(2.*velocity(nz-1,1).*delt + sqrt2*delx).*...
       (snapshot(nz-2,1)/delx + snapshot(nz-1,2)/delx +...
       sqrt2/(velocity(nz-1,1)*delt)*snap2(nz-1,1));

snapshot(nz,2) = velocity(nz,2).*delt.*delx./(2.*velocity(nz,2).*delt + sqrt2*delx).*...
       (snapshot(nz-1,2)/delx + snapshot(nz,3)/delx +...
       sqrt2/(velocity(nz,2)*delt)*snap2(nz,2));

snapshot(nz,1) = velocity(nz,1).*delt.*delx./(2.*velocity(nz,1).*delt + sqrt2*delx).*...
       (snapshot(nz-1,1)/delx + snapshot(nz,2)/delx +...
       sqrt2/(velocity(nz,1)*delt)*snap2(nz,1));

if boundary == 1

  % for upper right hand corner

  snapshot(2,nx) = velocity(2,nx).*delt.*delx./(2.*velocity(2,nx).*delt + sqrt2*delx).*...
         (snapshot(3,nx)/delx + snapshot(2,nx-1)/delx +...
         sqrt2/(velocity(2,nx)*delt)*snap2(2,nx));

  snapshot(1,nx-1) = velocity(1,nx-1).*delt.*delx./(2.*velocity(1,nx-1).*delt + sqrt2*delx).*...
         (snapshot(2,nx-1)/delx + snapshot(1,nx-2)/delx +...
         sqrt2/(velocity(1,nx-1)*delt)*snap2(1,nx-1));

  snapshot(1,nx) = velocity(1,nx).*delt.*delx./(2.*velocity(1,nx).*delt + sqrt2*delx).*...
         (snapshot(2,nx)/delx + snapshot(1,nx-1)/delx +...
         sqrt2/(velocity(1,nx)*delt)*snap2(1,nx));

  % for upper left hand corner

  snapshot(2,1) = velocity(2,1).*delt.*delx./(2.*velocity(2,1).*delt + sqrt2*delx).*...
         (snapshot(3,1)/delx + snapshot(2,2)/delx + ...
         sqrt2/(velocity(2,1)*delt)*snap2(2,1));
 
  snapshot(1,2) = velocity(1,2).*delt.*delx./(2.*velocity(1,2).*delt + sqrt2*delx).*...
         (snapshot(2,2)/delx + snapshot(1,3)/delx + ...
         sqrt2/(velocity(1,2)*delt)*snap2(1,2));

  snapshot(1,1) = velocity(1,1).*delt.*delx./(2.*velocity(1,1).*delt + sqrt2*delx).*...
         (snapshot(2,1)/delx + snapshot(1,2)/delx + ...
         sqrt2/(velocity(1,1)*delt)*snap2(1,1));

end